/* 
 * File:   soah_common.h
 * Author: boyan
 *
 * Created on April 18, 2008, 2:34 PM
 */

#ifndef _SOAH_COMMON_H
#define	_SOAH_COMMON_H

#include "glodef.h"
#include "vecs.h"
#include "fe_matr.h"

extern int dump_bin;
extern int MAX_PERIODS; // the maximum number of periods to compute periodic solution
extern int MAX_SHAPES; // the maximum number of steps of the shape optimization algorithm


#ifdef	__cplusplus
extern "C"
{
#endif


    void set_root_dir(const char *dir);
    const char *root_dir();
    const char *nodes_dir();
    const char *logs_dir();
    const char *state_dir();
    const char *adjoint_dir();

#ifdef	__cplusplus
}
#endif

extern void *solver_data;

extern vector eu, ev;
extern scalar ep, eq;


extern const unsigned elemU;
extern const unsigned elemP;

extern double alpha;
extern double beta;

extern double total_cost;


/* shape gradient and related quantities */

/* The following fields are defined only on the surface of the domain.  However,
 * for convenience, values for all nodes are allocated, but only values at
 * boundary nodes are significant 
 */
extern double sg_norm;
extern double sg_norm_inf;
extern double sg_lambda;
extern vector sg_z; /* the direction vector field */
extern scalar sg_xi; /* the size of diaplacement scalar field */
extern scalar sg_w1; /* auxiliary memory */
extern scalar sg_w2; /* auxiliary memory */
extern scalar sg_w3; /* auxiliary memory */

#define NORM_ERR_U_L2           0
#define NORM_ERR_P_L2           1
#define NORM_ERR_U_H1           2
#define NORM_ERR_P_H1           3
#define NORM_DUDT               4
#define NORM_U_L2               5
#define NORM_P_L2               6
#define NORM_GRAD_U_L2          7
#define NORM_EXCT_U_L2          8
#define NORM_EXCT_P_L2          9
#define NORM_EXCT_GRAD_U_L2     10
#define NORM_NUMBER             11

extern double current_solution_norm[NORM_NUMBER];
extern double time_inf_solution_norm[NORM_NUMBER];

extern double current_adj_solution_norm[NORM_NUMBER];
extern double time_inf_adj_solution_norm[NORM_NUMBER];

extern node_t *reference_nodes;
extern int *can_move_nodes;
extern fe_matr_p move_grid_matr;

#define USE_WEAK_GRADIENT   0
#define COST_FUNC_SQUARE    1
#define SHAPE_GRAD_SMOOTH   0

#ifdef	__cplusplus
extern "C"
{
#endif

    /* these are in soah_common.c */
    void set_defaults(void);
    void move_grid(void);
    double CostFunctional(void);
    void soah_check_heart_geometry(void);

    /* these are in soah_shape_gradient.c */
    void make_sg_z(void);
    void alloc_sg(void);
    void free_sg(void);
    int compute_shape_gradient(void);

    enum
    {
        sigma_open = FALSE, sigma_closed = TRUE
    };
    int _is_face_on_sigma(const grid3d_t *grid, int fn);
    int _is_node_on_sigma(const grid3d_t *grid, int nn, int closed);



    /* these are in soah_shape_optim.c */
    double soah_shape_cost(vector init_u, scalar init_p);
    void soah_shape_optimize(vector init_u, scalar init_p,
            vector finl_v, scalar finl_q);
    void soah_shape_optimize_armijo(vector init_u, scalar init_p,
            vector finl_v, scalar finl_q);
    void soah_shape_optimize_pit(vector init_u, scalar init_p,
            vector finl_v, scalar finl_q);

    void soah_shape_optimize_test_gradient_at_point(vector init_u, scalar init_p,
            vector finl_v, scalar finl_q);

    /* this is in soah_mesh.c */
    int soah_generate_meshes(int error);
    int soah_move_solid_wall(void);
    void set_reference_nodes(const node_t *new_nodes);

    /* these are in soah_state.c */
    void init_state_solver(void);
    int solve_state_ivp(vector u_init, scalar p_init)
    __attribute__((warn_unused_result));

    void done_state_solver(void);
    void print_conv_norms(void);


    /* this is in soah_adjoint.c */
    int solve_adjoint_periodic_fvp(const vector v_final, scalar q_final)
    __attribute__((warn_unused_result));
    int solve_adjoint_shape_fvp(const vector v_final, scalar q_final)
    __attribute__((warn_unused_result));
    int solve_adjoint_analytic_fvp(const vector v_finl, scalar q_finl,
            vector vtT_bc)__attribute__((warn_unused_result));
    void print_adj_conv_norms(void);


    /* this is in soah_periodic.c */
    /* ATTENTION: solve_state_periodic overwrites the input vectors */
    int solve_state_periodic(vector u_init, scalar p_init
            )__attribute__((warn_unused_result));
    int solve_shape_adjoint_periodic(vector v_finl, scalar q_finl
            )__attribute__((warn_unused_result));
    int solve_state_periodic_armijo(vector u_init, scalar p_init
            )__attribute__((warn_unused_result));
    int solve_state_periodic_adjoint(vector u_init, scalar p_init
            )__attribute__((warn_unused_result));
    int solve_state_periodic_lbfgs(vector u_init, scalar p_init
            )__attribute__((warn_unused_result));


    /* this is in soah_broyden.c */
    int solve_state_periodic_broyden(vector u_init, scalar p_init
            )__attribute__((warn_unused_result));
    int solve_adjoint_shape_periodic_broyden(vector v_finl, scalar q_finl
            )__attribute__((warn_unused_result));


#ifdef	__cplusplus
}
#endif

#endif	/* _SOAH_COMMON_H */

